#这个代码是为了对比ASM文件和WGS文件，如果WGS上该位点为杂合而且在ASM上找不到数据，就标记为FALSE
setwd("C:/Users/Administrator/Desktop/ASM/")
file1=read.table("interested_position_in_WGS2.txt",head=T,sep="\t",stringsAsFactors = F)
file2=read.table("interested_position_in_ASM.txt",head=T,sep="\t")
row.names(file2)=file2$unitID
header=c("M18","M20","M22","M40","M26","M27","M5","M36","M1","M30","M41","M47","M49","M51","M7","X1T","M43")

for(i in 1:nrow(file1)){
x=0
y=0
for(j in 12:13){
a=unlist(strsplit(as.character(file1[i,j]),":"))
if(a[1]=="0/1" |a[1]=="0|1"){
x=x+1
heter=paste(header[j-11],"_var_freq",sep="")
#unitID=paste("chr",file1[i,]$unitID,sep="")
unitID=file1[i,]$unitID
if(is.na(file2[unitID,heter])){
y=y+1
file1[i,j]=paste("FALSE",file1[i,j],sep="_")
}
}
}
file1[i,"All_heter"]=x
file1[i,"missing_heter"]=y
}



file1=read.table("E:/ASM/countdiff只考虑单发/interested_position/unique_id_61positions_in_M49-X2B_WGS.txt",head=T,sep="\t",stringsAsFactors = F)

row.names(file1)=file1$unitID
file2=read.table("61positions.txt",head=T,sep="\t")
row.names(file2)=file2$unitID
header=c("M49","X1T")
positions=as.character(file2$unitID)
file1=file1[positions,]


sel1=paste(unlist(strsplit(gsub(".snp","",sel),"_")),"_reads1",sep="")
sel2=paste(unlist(strsplit(gsub(".snp","",sel),"_")),"_reads2",sep="")